Inertial modes in stratified rotating neutron stars : An evolutionary description 
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With (nonbarotropic) equations of state valid even when the neutron, proton and electron content 
of neutron star cores is not in beta equilibrium, we study inertial and composition gravity modes 
of relativistic rotating neutron stars. We solve the relativistic Euler equations in the time domain 
with a three dimensional numerical code based on spectral methods, in the slow rotation, relativistic 
Cowling and anelastic approximations. Principally, after a short description of the gravity modes 
due to smooth composition gradients, we focus our analysis on the question of how the inertial modes 
are affected by nonbarotropicity of the nuclear matter. In our study, the deviation with respect to 
barotropicity results from the frozen composition of nonsuperfluid matter composed of neutrons, 
protons and electrons, when beta equilibrium is broken by millisecond oscillations. We show that 
already for moderately fast rotating stars the increasing coupling between polar and axial modes 
makes those two cases less different than for very slowly rotating stars. In addition, as we directly 
solve the Euler equations, without coupling only a few number of spherical harmonics, we always 
find, for the models that we use, a discrete spectrum for the Z = m = 2 inertial mode. Finally, 
we conclude that, for nonbarotropic stars, the frequency of this mode, which is our main focus, 
decreases in a non-negligible way, whereas the time dependence of the energy transfer between polar 
and axial modes is substantially different due to the existence of low-frequencies gravity modes. 



I. INTRODUCTION 



Since the result of Andersson and Friedman & 
Morsink numerous works have dealt with inertial 
modes of rotatingrelativistic stars (see reviews by Fried- 
man & Lockitch Q , Andersson Q and Kokkotas & RuofF 
[3). In perfect fluids, these modes have indeed been 
proven to satisfy the Chandrasekhar-Friedman-Schutz 
(CFS; and Q) criterion for instability whatever the 
angular velocity of the star is. Hence, due to their 
coupling to the gravitational field, these oscillations 
make neutron stars (NSs) very promising sources of 
gravitational waves (GW). Yet, as NSs are complex rela- 
tivistic objects, far from perfect fluids, the astrophysical 
relevance of the instability of inertial modes in rotating 
NSs is still an open issue. The final answer to this 
question will depend on the many physical phenomena 
that take place within NSs: their initial state at birth 
Ip, existence of differential rotation , of a crust 
(Ha, Hp, El, El), of a huge malefic field of 
superfluid states of nucleons (|l5|. flq . E3); exotic 
particles, etc. Understanding better that potential 
instability is important for two main reasons. First, 
if it were proven to be relevant, it could have a great 
impact on our idea about the evolution of NSs and might 
explain the relatively slow angular velocity of observed 
pulsars. Moreover, in order for the data collected by 
GW detectors like VIRGO or LIGO to be useful, they 
have to be compared with predicted signals, such as the 
signal that results from the instability of inertial modes. 
Since theoretical predictions depend on the models used 
for the inner structure of NSs, this confrontation could 



possibly be as instructive for nuclear physics at high 
density as was helioseismology for neutrino physics. 

In a previous article (T^ (referred to here as Paper 
I), we have presented a spectral hydrodynamical code, 
which uses spherical coordinates, with the aim of study- 
ing time evolution of inertial modes in slowly rotating 
neutron stars. Having discussed the motivation for the 
slow rotation approximation in this previous article (it is 
still good even for the fastest known pulsars, see Paper 
I), we keep working with it in the present paper. Thus, 
the next step in our project was to improve the linear 
study of relativistic inertial modes by trying to take into 
account in a more realistic way the microphysical con- 
ditions that occur inside actual NSs. Among them, we 
retained the stratiflcation of the star. It was quite easy 
to implement, and it is well-known that inertial modes 
change drastically depending on whether the equation 
of state (EOS) that describes the NS is barotropic or 
nonbarotropic^, while the status of purely axial modes 
of relativistic nonbarotropic NSs is still under debate 
([3,I3)- In addition, it is easy to show that, for the 
most basic model of a star built with the so-called "npe" 
matter, the assumption of a barotropic EOS, even if 
easier to use, is wrong as soon as the star is cold enough 
and the matter not in a stationary state. Hence, here we 
consider the time evolution of inertial modes in slowly 



^ We remind that an equation of state is said "barotropic" if the 
pressure P depends only on one variable, for instance the bary- 
onic number density nj. 
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rotating relativistic NSs composed of neutrons, protons 
and electrons, with an "out-of-equilibrium" EOS. It 
should be kept in mind that in the following, we shall call 
"nonbarotropic" a peculiar type of nonbarotropic EOSs: 
the situation of npe matter with frozen composition. 

The plan of the paper is as follows. In Section we 
start with an overview of the basic assumptions made in 
the study of the structure of cold neutron stars, with em- 
phasis on the beta equilibrium between neutrons, protons 
and electrons, which constitute, for not too massive NSs, 
the most important part for studies of their oscillations: 
their core. We end this section with a description of our 
calculation of stationary rotating NSs and of the EOS 
used to get configurations to perturb in the following sec- 
tions: the Prakash-Ainsworth-Lattimer (PAL) EOS [l9l |. 
In Section IIIII a short summary of the microphysics of 
"perturbed" NSs composed of npe matter is given. In this 
section, we justify the assumption of frozen composition, 
which is then used in Section Hvl that describes the basis 
of the linear hydrodynamical study. The main feature of 
a nonbarotropic star is the existence of a nondegenerate 
spectrum for its polar modes, even when it is not rotat- 
ing. As frozen composition results in a nonbarotropic 
effective EOS, these so-called gravity modes ((7-modes) 
are studied for the PAL EOS in Section ^ After which, 
in Section lVll we describe inertial modes of rotating NSs 
with the PAL EOS but without the (correct) hypothe- 
sis of frozen composition, before ending with gravity and 
inertial modes in rotating NSs taking into account the 
frozen composition. Finally, the results are discussed in 
the Conclusion. 



II. FULLY CATALYZED ROTATING 
RELATIVISTIC NEUTRON STARS 

A. Model of the neutron star matter 

The outer layers of cold NSs are known to play key 
roles in several physical phenomena, as the early and 
late cooling of pulsars |23|, the existence of glitches 
((U, [23, ^) and the spectra of X-ray bursts 
Yet, since more than 95% of the mass of cold NSs is 
included in their core, outer layers are not relevant at 
the most basic levels of studies of both the structure of 
rotating relativistic stars and their oscillations. Thus, in 
the following, a NS is modeled by its core, and the outer 
layers are only discussed as time independent boundary 
conditions for the oscillations (see Section . 

The minimal nontrivial description of a core of NS 
is a mixture of neutron, protons, and electrons at a 
temperature that is sufficiently low for the matter to be 
degenerate but sufficiently high for the nucleons not to 
be superfluid. Such conditions, which imply that thermal 
contributions to the EOS can be neglected, are satisfied 
for 10^ K <T < 10^" K, at densities above one tenth of 



normal nuclear density no = 0.16 fm~ . Moreover, when 
a NS has cooled down enough for these assumptions to 
be valid, the matter has become transparent for neutri- 
nos, which therefore do not contribute to the EOS either. 

Another physical constraint to take into account for 
the study of both unperturbed and perturbed stars is 
the electrical charge neutrality, which comes from the 
very short characteristic time scale of electromagnetic 
processes. In degenerate npe matter, the so-called plasma 
frequency of electrons, whose inverse give the typical time 
for the plasma to maintain charge neutrality, is 



iire'^neC^ 

Me 



,2\ 1/2 



(1) 



where fie is the chemical potential of electrons similar to 
fie ^ 122 (ne/no)^/^ MeV (see below), with the electron 
density rig. This leads to an electromagnetic relaxation 
time Tpe = 4 X 10~^^(ne/no)^^^ s, which implies that, 
for our purpose (phenomena with typical time scales of 
some ms), matter can always be considered as neutral, 
and there is equality between the electron and proton 



densities 



This relation also makes thermody- 



namical quantities functions of nucleon (baryon) density 
nb = nn + np and proton fraction Xp = np/n^ only. Fur- 
thermore, at baryon densities of the order of ng, nucle- 
ons form a strongly interacting Fermi system and matter 
is very rich in neutrons, while we can replace nucleon 
masses (which are assumed to be equal) by the neu- 
tron mass = m„ = 939.57 MeV. If we neglect now 
the Coulomb interactions between charged constituents 
of matter and treat the electrons as a free ultrarelativistic 
Fermi gas, the energy per nucleon (excluding rest energy 
of nucleons) can be written 



E{nb,Xp) = E^{nb,Xp) + Ee{nb,Xp) 



(2) 



where Ejq is the nucleon contribution, and the electron 
term Ee is 

3 1 i 

Ee{nb,Xp) = -bu^x^ , b = hc{^TT'^noY/^ = 331.4MeV, 

(3) 

with u = nb/n^. 

Extensive many-body calculations with realistic two- 
nucleon and three-nucleon forces show that in the whole 
range of < a;^ < 1, the energy per nucleon is (to a very 
good approximation) quadratic in the neutron excess pa- 
rameter (n„—np)/n;, = 1 — 2a;p ([23, l^^). At ~ 0.5, 
characteristic value for terrestrial nuclear physics, this is 
just the consequence of the charge symmetry of strong- 
interaction. At a;p ~ 1, this is a numerical fact. There- 
fore, the a;p-dependence factors out, and E^ is given by 
the formula 



E^{nb,Xp) = Wo{nb) -l- S{nb){l - 2xpy 



(4) 



Here, Wo('^h) is the energy per nucleon of the symmet- 
ric nuclear matter {xp = and S{nb) is the symmetry 
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energy at nucleon density rifc. Many possibilities exist 
for these two functions, among them those that form the 
PAL EOS This EOS, which is retained in this study, 
will be described in more detail with the presentation of 
the obtained background stars used for our modes cal- 
culations f Section III D(l . However, to be able to calcu- 
late these unperturbed stars, another relation is needed 
to close the system of equations. This relation comes 
from the hypothesis of beta equilibrium between neu- 
trons, protons and electrons. 



B. EOS in beta equilibrium 

For matter transparent to neutrinos, the equilibrium 
with respect to the weak-interaction beta-processes 



P 



p + e 



n + l>e 



(5) 



implies a relation^ between the chemical potentials of the 
matter constituents, fij ^dp/drij, with 



p{nb,Xp) = nb{mNC^ + E) 



(6) 



the total energy density [energy per baryon E that comes 
from Eq.|(21) plus nucleon rest energies], and j = n,p,e. 
For matter transparent to neutrinos, this relation reads 



fj.p + fie 



(7) 



With our quadratic form of the a:p-dependence of -En, 
Eq.Q), the beta equilibrium condition involves only the 
nuclear symmetry term. We have 



fj.p- fie = 4 5'(1 - 2Xp) 



bu^x^ . (8) 



The resulting fraction of protons for matter in beta 
equilibrium, denoted by x'^{nb), is determined from 



ay 



45' 



(9) 



(10) 



where y^^xl'i and a{rn,) = [bu^/^) / S). 

The last equation admits analytical solutions (c/. Car- 
dan's method) among which the only real one is 



1 _ oM 
y^b) 2 22/3 



(FK) -1-1)1/3 _ (FK)- 1)1/3 

(11) 



with T{nb) - 



1 
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Then, the barotropic EOS in beta equilibrium, which is 
used for the calculation of unperturbed stars, is obtained 
from 



Pcq{nb) = p{nb,x°^{7n,)) 
with the total pressure 



2 9i Pcq/nb) 
dub 



(12) 



(13) 



Note that when the chemical potential of the electrons 
is above the rest mass of muons, rn^ ^ 105.65 MeV, tak- 
ing into account the appearance of the latter via the equi- 
librium between the processes 



and p 



(14) 



would have implied higher proton fractions. But as a 
rule, the role of muons is practically irrelevant, as far as 
the EOS is concerned, and we shall always neglect them. 



C. Stationary rotating relativistic neutron stars 

The numerical solution of the Einstein's equations for a 
cold stationary rotating relativistic star with a barotropic 
EOS, P{p) [reached by Eq.lO and (HSJ], has now be- 
come a classical problem (see for a review) whose 
detail we will not enter into here. Let us just remind the 
reader of the main assumptions: 

- there are two Killing vectors k (timelike at spacelike 

infinity) for stationarity and m (with closed orbits, 
spacelike everywhere except on the so-called rota- 
tion axis where it vanishes) for axisymmetry; 

- spacetime is asymptotically flat, the metric g "tends 

toward" the flat Minkowski metric r\ at infinity; 

- matter is a perfect fluid whose energy-momentum ten- 

=3 



sor is'^ 



T=(p -t- P)u( 



(15) 



where u is the 4-velocity of the fluid. Obviously, 
P and p are measured in the frame comoving with 
the fluid and thus are the equilibrium functions ob- 
tained in the previous section; 

the velocity field corresponds to rigid axisymmetric ro- 
tation and can be written 



u = A(k -I- Jim) 



(16) 



Note that this relation remains the same even when direct Urea 
processes, presented here, are not allowed due to kinematical rea- 
sons, and have to be replaced by processes that involve spectator ^ Except where otherwise noted we use natural units in which 
nucleons, the modified Urea reactions. c = Qj^ = 1. 
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where fl is by definition the angular velocity of the 
star and A a scalar field such as the norm of u is 
— 1 (with our convention for the signature) . 



With these assumptions, the energy 
satisfies the condition of circularity [25 



momentum tensor 
I and the maximal 



slicing quasi-isotropic coordinates [2£ 
write the metric as 



can be adopted to 



{N^ - N"^) dt^ ~2N^dtdip + 

I 



A* 



(dr 



dif + B'^ sin^ d di/ 



(17) 



where we have introduced the notation of the 3+1 formal- 
ism: N is the lapse and the third component of the 



shift 3- vector, with 



^4 _g2 ^v^2 gjj^2 ^ 



the covariant (^-component of the latter. 

Because of the existence of two Killing vectors all the 
functions only depend on (r , i9), and this choice of coor- 
dinates reduces the Einstein's equations to a system of 
elliptic equations that we solve with the fully relativistic 
spectral code presented in Bonazzola et al. f^. This 
code also provides us with what will be the velocity field 
to perturb in the study of oscillations (Section II V|). In- 
deed the A scalar field, which normalizes the velocity, can 
be shown to satisfy 



A 



(|k-k + 2r2k-m + f^^m-ml)" 



(18) 



where • denotes the scalar product associated with the 
metric. 



D. Calculated background stars 

In order to study inertial modes of NSs more realistic 
than simple barotrops, we need an EOS for the nuclear 



J 



(0), 



2^" + - 



Bu" 



B' u" 



matter that is also valid when the fluid is not in beta 
equilibrium (O . Indeed, as explained in the Introduction 
and as calculations will show in Section IIIII a lump of 
matter transported out of its equilibrium position by 
an oscillation will not adjust its composition to its new 
environment instantaneously. This will require some 
time, which has to be taken into account in the dynamics 
since it allows new modes (see Section Hvj) to exist. 



Among the numerous EOSs that describe nuclear mat- 
ter in NSs, we decided to use the equation proposed by 
Prakash et al. [Toj which is among the few which are 
valid also out of beta equilibrium. Moreover, this analyt- 
ical EOS based partly on the experimental properties of 
nuclear matter, and partly on the results of many-body 
calculations of asymmetric nucleon matter, is in fact a 
set of EOSs. Indeed, several possibilities are given for 
the choice of the symmetry energy function Sijih), while 
the energy per nucleon of the symmetric nuclear matter 
Wo(?T.fc) depends on some parameters. Thus, following 
Prakash et al., we have 



E 

i =1,2 



PF 

A. 



—— — arctan 



EL 
A. 



(19) 



where 



E 



is the free nucleon gas Fermi energy at saturation 
density (~ 36 MeV); 

the dimensionless variable u is u^nt / no, where 
no = 0.16 fm^"^ is the saturation density of sym- 
metric nuclear matter; 

the second and third terms, that depend on the con- 
stants A, B, B' and cr, reproduce the static inter- 



actions between nucleons and are chosen in such a 
way that the EOS is causal; 

the last two terms, that depend on the momentum, 
mimic the dynamical part of the strong interaction. 
In those terms, the A parameters are characteristic 
lengths associated to coupling constants Ci. The 
first term (such as Ci < and Ai = 1.5pF°) is 
attractive at large distances, whereas the second 
(C2 > with A2 = 3^^-") reproduces the re- 
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pulsive behavior of the strong interaction at short 
scales. 

In order to separate the kinetic and potential contri- 
butions in the symmetry energy, Prakash et al. 19} write 

5K) = (2i - l) ^4"^ (ui - F{u)) + SoF{u), 

(20) 

where F(u) is a function which describes the 
(poorly known) density dependence of S. To get 
S{no) = 5o = 30 MeV, one has F(l) = 1. 

Playing with the form of the functions Wo{ni,) and 
S{nb) enables one to change the values of the com- 
pression modulus and the symmetry energy. However, 
this analytical EOS, that fits known results around 
the saturation density, has the drawback that, at low 
densities, it is no longer physical. Indeed, to be realistic 
at low densities, any EOS has to take into account the 
appearance of nuclei. Hence, for several choices of the 
key quantities (such as the value of the compression 
modulus and the symmetry energy) within the proposal 
of Prakash et al., pathological features like negative 
pressure or convective instability appear. Moreover, 
even for the EOSs that do not present these awkward 
characteristics, the low density part is unphysical and 
some prescription has to be adopted for the study of 
modes. 

Among the functions and the set of parameters that 
Prakash et al. propose, we retain F(u) = \[u and take 
values of the parameters that give a compression modu- 
lus of either 180 MeV (model A) or 240 MeV (model B). 
In addition, to prevent the troubles linked to low densi- 
ties, our choice is to cut the star below a given threshold 
density. As discussed in Section ^1 since the mass of 
the crust is quite small, this cutoff does not really af- 
fect either the metric or the main properties of the star. 
Moreover, adding some boundary conditions on the sur- 
face described by the value of the density chosen for the 
cutoff is an easy way to mimic the phase transition that 
occurs between the outer-core and the inner-crust. This 
point will be discussed in more detail in the sections that 
deal with modes (Sections |E1 El ED but it was noticed 
now to clarify Tabled Indeed, in this table are summa- 
rized the background configurations that we use. But in 
order to calibrate the effect of the cutoff on the modes, 
we tried several values for the density threshold. Hence, 
in Tabled we show structural properties of stars for the 
two EOSs that we keep with different central densities 
and the values that we take for the cutoff density**: 0.08 
or 0.1 baryon.fm"'^. As a first validation of the cutoff 



procedure, we can check that for given EOS and central 
density, the cutoffs imply differences in the gravitational 
mass of some few percent in the worst case. Finally, note 
also that the central (maximal) proton fraction is always 
below 1/9, which forbids direct Urea processes for npe 
matter (see also Section llll C|) . 



III. PERTURBED CONFIGURATIONS I: 
MICROPHYSICS 

A. npe matter 

When trying to deal with modes of NSs in a not too 
unrealistic way, i. e. when trying to go beyond the simple 
model of a single barotropic fluid, the most basic model 
to consider is a NS composed of neutrons, protons and 
electrons. Remember that we deal here only with NSs 
sufficiently old to have become completely transparent 
to neutrinos, but we neglect possible superfluidity of 
nucleons. Anyway, even without neutrinos to complicate 
the game, these three fluids could, a priori, be only 
partially coupled and then three relativistic Euler 
equations (with coupling terms) could be needed^. 

But, as discussed in Section^ due to the very short 
time scale of the Coulomb interaction and due to its long 
range, it is well-known that, for oscillations with periods 
typical of NSs modes (some ms), protons and electrons 
are in fact interdependent with the same velocities and 
densities: matter is locally electrically neutral Xp = Xe- 
The first thing that we shall discuss now, is how the 
strong interaction between nucleons plays a similar but 
not exactly identical role. It makes the velocities of pro- 
tons and neutrons equal (without superfluidity), while 
not leveling out their densities, as the latter depend on 
the weak- interaction (see Section IIIB|I . As the demon- 
stration of this in the most general case is quite tricky, we 
shall work here with a simpler situation that provides an 
order of magnitude estimate sufficient for our purpose. 



B. Damping of the relative n — p flow 

We assume that there is no external force acting 
differently on n and p (i.e., we neglect the effects of 
the electromagnetic field). Consider the relative motion 
of neutrons and protons along the z-axis, with no 
baryon number flow, Vn = —Vp = v. Neglect all spatial 
gradients, and assume, that at t = we have w(0) = vq. 
By solving the kinetic equations for neutron and proton 



* i.e. one half or almost two thirds of the saturation density no- 
This range is consistent wi th r ealist ic e valuations of the density 
at crust-core interface (see l3Cl and ISlI ') 



^ Notice that in this case, the easiest way to reach the equations of 
motion would not be to start with the usual conservation of the 
energy-momentum tensor, but to us e th e Lagrangian formalism 
for multifluids developed by Carter |32| . 
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TABLE I: Properties of the NSs that are used as backgrounds for the calculation of inertial and gravity modes. The entries in 
the table are: type of EOS (name of model); central relativistic enthalpy (dimensionless) ; central baryonic density (in fm~^); 
central proton fraction; density at the surface; enthalpy at the surface; gravitational mass (in solar mass units), circumferential 
equatorial radius (in km) and compactness C = Ma/Reg (in dimensionless units). Note that Model B is stiffer and allows a 
minimal gravitational mass which is larger than 1.2 Mq. This is the reason why we do not have stars of model A and B with 
the same 1.2 Mq mass, whereas we have stars of the same compactness A(26) and B(24), but also of the same 1.6 Mq mass 
A(4) and B(3). 



EOS 


Model A(21) 


Model A(26) 


Model A(4) 


Model B(24) 


Model B(3) 


he 


0.21 


0.26 


0.4 


0.24 


0.3 


ric [fm^^] 


0.573 


0.676 


0.966 


0.540 


0.635 


^Pc [%] 


7.43 


8.05 


9.47 


7.22 


7.81 


nsurf [fm~^] 


0.08 0.1 


0.08 0.1 


0.08 0.1 


0.08 0.1 


0.08 0.1 


/lBurf[10-"] 


2.0 2.3 


2.0 2.3 


2.0 2.3 


1.9 2.2 


1.9 2.2 


Mg [Mq] 


1.20 1.18 1.17 


1.35 1.34 1.33 


1.60 1.59 1.59 


1.42 1.41 1.40 


1.60 1.59 1.58 


Req [km] 


11.3 10.4 10.3 


11.1 10.4 10.3 


10.2 10.0 9.9 


11.8 11.1 11.0 


11.6 11.0 10.9 


C [10-1] 


1.57 1.67 1.68 


1.79 1.91 1.91 


2.33 2.36 2.38 


1.77 1.88 1.89 


2.03 2.13 2.14 



distribution functions, with a collision term due to the 
strong-interaction nucleon-nucleon scattering, we will 
find that v{t) = z)(0)e-*/'^"'' , where r„p is the relaxation 
time after which the relative n — p flow is damped. As 
a nuclcon fluid forms a strongly interacting dense Fermi 
liquid, the n — p scattering cross section to be used in 
the kinetic equation is different from that in vacuum, 
and the transport process refers to the quasiparticles of 
the Fermi liquid theory (see, e.g., Baym et al. [2l|). 

For symmetric nuclear matter at normal nuclear 
density, one gets, using the Fermi liquid theory, 
Tnp ~ 10-1^ Tg-^ s (Haensel where Tg is the tem- 

perature in 10^ K. For a highly asymmetric neutron star 
matter, with Xp = 0.05, simple estimate in which the 
vacuum n — p scattering cross sections are used, were ob- 
tained by Baym et al. |23, T„p ~ 10"^^ Tg"^ s. Hence, in 
view of (microscopic) smallness of r„p compared to any 
other time scales relevant for our study, neutrons and 
protons can be assumed to move together as one single 
nucleon fluid. Yet, as already mentioned, the ratio be- 
tween their densities is not governed by strong interaction 
but by weak-interaction, which makes their relation dif- 
ferent from the relation between electrons and protons as 
we shall see now. 



C. Nonequilibrium beta-processes 

Since the time for beta equilibrium ^ to settle in de- 
pends on the production of thermal neutrinos, for the 
time being we shall consider an element of hot npe mat- 
ter at temperature T . This lump has an instantaneous 
nucleon density = rip -!-«„, which due to its temporal 
evolution (that results from the matter flow) is out of 
beta equilibrium, so that 

= - M« 7^ 0. (21) 



Let us also characterize this deviation from beta equilib- 
rium by ^ = Xp — Xp^{ni,). 

Each of the three fluids (neutrons, protons and elec- 
trons), is separately in thermodynamic equilibrium, due 
to strong and electromagnetic interactions. However, 
weak-interaction processes may be too slow to establish 
beta equilibrium at [nf,,T]. We denote the values of n„ 
and Hp, which correspond to beta equilibrium at [nti,T], 
by n';^, n'p. The values n^, n'^ are the reference ones, 
for a given [n(,,T] pair, and are compared with actual 
values n„, Up. 

A nonzero value of (5/i, Ea. (|21|l . implies nonequilib- 
rium reactions, which tend to decrease deviation from 
beta equilibrium (in accordance with the Le Chatelier 
principle) . Let us denote the rate of the change of Up at 
flxed Hh by AT{nb,T,6fi) (see Sect. 3.5 of Yakovlev et 
al. We then have 

AT{nb,T,Sfi = 0) = , 
Ar(nb, T, -6(1) = -Ar(nb, T, S(i) , (22) 

while at flxed rif,, we can write 

rip = -n„ ^ = Ar(n(,,r, 5/z) , (23) 

where the dots denote time derivatives. 

In the linear approximation, ^ <^ x'p, 5^ implied by f 
at a flxed is 

where all cocfflcients are calculated at nii,Xp'^. 

In what follows, we assume in addition 
that Sfi/{T:kT) <^ 1, which is equivalent to 
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<C 0.27 Tg. This means that matter de- 
viates only weakly from the beta equilibrium, an 
approximation which is valid for small amplitude of 
neutron star pulsation. Moreover, we remind the 
reader that nucleons are supposed to be normal (no 
superfluidity, valid for Tg > 1) and we will consider 
here only cases with Xp < 1/9, which is in perfect 
agreement with our background models (see Table QJ, 
so that only modified Urea processes are allowed (see 
Lattimer et al. HJl). Notice that the case of direct Urea 
processes would anyway lead to a similar conclusion, just 
changing quantitatively the result. Under the previous 
assumptions, the formula for AT reads, keeping only 
terms linear in 



where normal nuclear density rig — 0.16 fm 



(25) 



In the relaxation time approximation, we can rewrite 
the formula for AT as 



Ar = 



rip - n'p'i 



T0 



(26) 



At fixed TT-b, relaxation from a state initially off equi- 
librium by 6n^ = — n'^, towards beta equilibrium, 
proceeds according to 



(27) 



Using the relation between 5^ and Sup, one shows that 
the beta-relaxation time is given by 



TP = 5.1 X 10^ (Tg)-^ — 



6 / '^h 



no 0.01 



1/3 



(28) 



Under the prevailing conditions, this time scale is so 
long compared to other typical time scales that it is a very 
good assumption to consider that every lump of matter 
keeps a frozen composition (value of Xp) in perturbed 
configurations^. In Section IIVI we shall see how, from 
the hydrodynamical point of view, this condition can be 
taken into account and leads to a nonbarotropic equation 
of state. Before that, we shall briefly discuss for npe 
matter viscosity coefficients and viscous damping times 
of oscillations. 



^ If hyperons were included, they would also move together with 
the bulk matter due to the strong interaction, but concerning 
the respective ratios, the situation would be a little different 
due to nonleptonic strangeness violating reactions. These weak- 
interaction reactions without Icptons like 



NA ^ Nn, nn ^ pS" 



(29) 



where N is a nucleon, have relaxation times of the order of the 



millisecond for T ~ 10^ K (t 



I). Thus the composition 



D. Viscosity of npe matter 

As we shall see later (Sections |3 and IVI|I . we deal 
in this study only with rotational and gravity m = 2 
modes, with the idea of focusing on the gravitational 
driven instability of the r-modes. For them, it is well- 
known 5] that shear viscosity kills the CFS instability 
(see Section IVI|) for low temperatures, whereas bulk- 
viscosity does it for high temperatures. On the other 
hand, gravity modes of cold NS star are not directly 
relevant for gravitational waves emission due to their 
very weak coupling with the gravitational field. Hence 
our goal here is just to "verify" that to use the equations 
of motion for perfect fluids is quite reasonable for the 
durations of our simulations (of the order of 1 second) 
during which viscosity would not really have time to 
damp the modes. 

The dissipation connected with viscosity of the npe 
matter is characterized by two density and tempera- 
ture dependent parameters, shear viscosity 77 and bulk- 
viscosity Shear viscosity results from the momen- 
tum transfer in the scattering processes between the npe 
matter constituents. It is mostly determined by large- 
momentum transfer collisions, and therefore the domi- 
nant contribution to rj is that of nonsuperfluid neutrons. 
An approximate analytic expression describing results 
obtained by Flowers & Itoh 37] is 



1.6 X 10i6(pi4)2 (rg)-2 g cm-i s-i 



(30) 



where pi4 = p/lO^'* g cm ^. 



On the other hand, bulk-viscosity results from the 
nonequilibrium beta-processes previously discussed. But 
remember that for the notion of bulk- viscosity to be valid, 
the deviations from the beta equilibrium measured by 
dfi should be much smaller than k T . In such condi- 
tions, by considering periodic pulsations of local pres- 
sure in the npe matter, one can obtain an expression for 
the mean heat deposition due to the nonequilibrium pro- 
cesses. When the relation between Fermi momenta 



PFn > PFp + PFe 



(31) 



is fulfilled, which in the case of npe matter corresponds 
to Xp < 1/9, the direct Urea processes are prohibited. 
As already mentioned above with the description of our 
background stars, this is always the case in this study 
(see Table U). Hence here only modified Urea (mUrca) 
processes are allowed, with their neutron branch 



n + n - 
n+p 



p + e + Ve , 
> n A- n + , 



and their proton branch 



could be no longer frozen for millisecond oscillations. 



p + n — > p + p + e + 
P + P + e > p + n + Vf. 



(32) 



(33) 
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The formulae obtained by Haensel et al. (see |3a|.|3g| 
for more detail) and applied to the npe matter give for 
the resulting bulk-viscosity 



mUrca 



lO^V'lrg)' g cm-1 s-i . (34) 



Here, lo/^ is the angular frequency of pulsations mea- 
sured in the units of 10"* s^^. Even if wc shall not 
need it here, let us mention that the magnitude of the 
bulk- viscosity increases by many orders if the direct Urea 
(dUrca) processes are allowed. These processes lead to 



/^dUrca 



1025^4 2 (r9)4e„pe g cm-1 s-1 , (35) 



where the threshold factor Qnpe vanishes if 

PFn > PFp + PFe- 

With Ea. H3U|) . we obtain the time scale for the shear 
viscosity j40j 



ode ^ X 10^ s , (36) 
Pl5 



where R = Rq x 10^ cm is the core radius, 
P — Pi5 X 10^^ g cm~^ the "mean" core density, 
and a^ojp a numerical coefficient that depends on the 
pulsational mode. 

Eq. I|35|l give a very rough estimate of the bulk- viscosity 
damping time scale due to mUrca bulk-viscosity: 



mU 



• i^mV ~ "mode J* 



(37) 



We see with these formulae that around 10^ K, none 
of the viscosities damps the modes very quickly^. More 
precisely, for Tg = 5 the damping time scale is of the 
order of hours. Hence for an evolution lasting a few 
seconds the assumption of perfect fluid is not so bad. 

As a final remark, let us mention that, in a star of 
radius R, the presence of a core of radius Rd where the 
direct Urea process is allowed would imply a damping 
time scale given by 



pR^ 

*modc ^dU 



C.dUrca ^ -, „3 
"mode ^ 



(38) 
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with, for the r-modes, a''''^^ much greater than one (see 
|4ll | for the case of the CPS instability). Additional in- 
crease of a^ad^'^ results from the fact that the sphere 
in which energy is dissipated has Rd for radius, that is 
only a fraction of the stellar radius R in which the mode 
propagate. However, the direct Urea process would not 
be for the modes as dramatic as the presence of hyperons 
and at Tg = 5 we would still expect damping on a time 
scale of minutes. 



the expected value of for the r-mode is significantly 

greater than one. 



IV. PERTURBED CONFIGURATIONS II: 
HYDRODYNAMICS 

A. Assumptions for the linear hydrodynamical 
study 

Doing time evolutions in general relativity is much 
more complicated than the solution of Einstein's equa- 
tions for stationary and axisymmetric configurations. 
Hence, it is natural (at least in a first study) to begin 
with some approximations in order to simplify the prob- 
lem. Thus, for the rest of this work, we assume that 

- the perturbed star is "not too far" (see for a proper 

definition of this) from the unperturbed one. It is 
then legitimate to linearize all equations (that gov- 
ern both the material fields and the gravitational 
field) with respect to a parameter that indexes a 
continuous family of time dependent solutions and 
whose value is for the unperturbed configuration. 
Obviously, this approximation means that in the 
following study of modes, the zero order terms are 
input that come from the unperturbed calculation 
presented in Section FH DI 

- the Eulerian perturbations of the gravitational field 

(i.e. of the metric) are quite "small" and not fun- 
damental for the physics. We then neglect them 
and write 



(39) 



where Sgfj^i, are the Eulerian perturbations of the 
metric (/^i,. This is the so-called relativistic Cowl- 
ing 43 approximation, introduced by Mc Dermott 
et al. 43] in the study of oscillations of warm neu- 
tron stars. As discussed by several authors, this 
approximation has several drawbacks in the rela- 
tivistic case ([3jE3|)j but it was also shown that it 
does not change dramatically the main feature of 
inertial modes (see |^ for a review). In addition, 
implementing the evolution of some metric pertur- 
bations (see Section lyij) in our code would require a 
non-negligible amount of work since boundary con- 
ditions have to be treated carefully with spectral 
methods 45] ; 

the NS is slowly rotating. Only terms linear in the 
angular velocity fl are kept in the equations, and 
every zero order term is assumed to only depend 
on the radial coordinate. This assumption was dis- 
cussed in Paper I, and we shall stress here just one 
point : its use makes it sufficient to work with Har- 
tle's equations 46] and not with the full Einstein's 
equations for the unperturbed calculation. Yet, as 
we plan to get rid of it in the future and since we 
have a fully relativistic code, it was easier for us to 
use this code to get the zero order terms. More- 
over, we verified in this way the relevance of the 
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slow rotation approximation. For more detail, see 
Paper I; 

- the spatial part of the metric is conformally flat 
(Ise nberg-Wilson-Mathews approximation, ^3 1 
|48| and [43), i.e. the metric ((T7|) becomes 

ds^ ^ -{N^ - N"^) dt^ - 2 TV,^ dt d^ + h^J dx' dx^ 

(40) 

with i, j e (1,2,3) and 

h = ar], (41) 
where a is the conformal factor and t] the flat Eu- 



where any function N, N'^ and a only depends on the ra- 
dial coordinate x^ = r, while dl^ is the length interval in 
the flat 3-space. With this metric, the linearized equa- 
tions of motion can be easily written in a way similar to 
the Newtonian equations. 

B. Equations of motion 

For reasons described in Section ITTll we assume that 
there is only one fluid whose relativistic equations of mo- 
tion are obtained by the usual conservation of the energy- 
momentum tensor of a perfect fluid 

V • T = , (43) 

with 

T = (p + P)u® u + Pg, (44) 

where V is the covariant derivative associated with the 
Levi-Civita afhne connection for the given metric H42|) . 



uf'it, r,'d,(j) 




Here, 5U^ , SU"^ , SU"^ and SV^ are first order quanti- 
ties added to the unperturbed 4-velocity that describes 



clidean 3-metric. The conformal approximation 
was shown to be very good for (even fast) rotating 
isolated NSs by Cook et al. |^, and once again 
using the fully relativistic code we could verify this 
before making the assumption in the linear part of 
the work. Moreover, as we already use the Cowling 
approximation, we do not really care to "kill again" 
the GW content of spacetime, which can at a first 
level be obtained with some "post-Newtonian mul- 
tipolar scheme" (see for instance [^). 

With all these assumptions, the final frozen metric 
used in the first order calculations is written 



(42) 

I 

For the isentropic motion of a perfect fluid with a given 
EOS, it is well-known that, due to thermodynamics, the 
system of equations formed by the baryonic number con- 
servation and the energy-momentum conservation l|43|) 
is degenerate. Hence, we shall work with the baryonic 
number conservation, 

V • (nfc u) = , (45) 

where rif, is the baryonic density measured in the fluid's 
local rest frame, plus the projections on the 3-space 
of the relativistic Euler equations (EE) obtained from 
Ea. (|I5|l . More precisely, the useful equations are those 
reached after linearizing close to a solution that describes 
the rigidly rotating star, Eqs. H16I18|I . 

As was briefly described in Paper I, in order to get 
equations as similar as possible to the Newtonian EE, 
we use the well-known results about rigid rotation of rel- 
ativistic stars, Eqs. Hlt)ll8|) . and write the perturbed 
4-velocity as 



(46) 



I 

rigid rotation, with 5U'^ that is not a dynamical variable 
but determined according to the constraint that u is a 



ds^= -[N'^ - r'^ sm[^f Nf"^) di^ ~ 2 sin[i?]^ iV^dt d(^ + dl^ 



l-h(5C/°[<,r,i?,0] 

{N[r] / a[r]f SVlt, r, t9, 0] 

{N[r] / a[r]f 5U'^[t,r,d,(j3\ / r 

VL + {N[r] / a[r]f W'P[t,r,d,(l)] / {rsm[d]) 
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4- velocity : u-u = — 1 for the chosen signature of the 
metric. Furthermore, SU^,5U''^ and JC/'^ are not con- 
travariant components of a 4- vector but convenient vari- 
ables used in our calculations, which are nevertheless the 



components of W, a 3-vector, in the orthonormal basis 
associated with the spherical system of coordinates for 
the flat 3-space. With all this, we obtain 



{dt + n 84.) SW + S [drP I f] - 2 SV^ sin[i?] {{n - N^) {r a' / a + 1 - r N' / N) - r N'^' / 2^ ^ , 



{df + n 84,) SU^ + S [diiP I (r /)] - 2 JC/'^ cos[i?] (f^ - iV^) = , 



(47) 



(dt + ^ d^) + 8 [d^P /{fr sin[i?])] + 2 sin[i?] (^{n ~ m) {r a' / a + 1 - r N' / N) - r Nf' / 2) 

+ 2 5U'^ cos[z9] {n - 
V 



where the ' indicates the derivation with respect to the 
radial coordinate, S is the operator for the Eulcrian 
perturbation and f = p + P. 

To this system of equations, we add the linearized 
baryon number conservation that, in the slow rotation 
approximation, can be written as 



{dt + na^) 5n + 



j\t2 , . 

-2- div \n\^^ , (48) 



where we define h = n\,N\r^ a\r\ and introduce the 
already mentioned 3-vector 14^, using also the usual 
3-dimensional spatial divergence operator. 

Hence, with the introduction of the relativistic anelas- 
tic approximation defined in Paper I (see also this paper 
for a discussion of the motivations for this approxima- 
tion), we finally get the Newtonian- like equation 



div 



0, 



(49) 



which is the equation used in the current version of the 
code. The only missing element to close the system of 
equations is now the way thermodynamical quantities 
(mainly the pressure) are perturbed, i.e. we need to give 
a prescription to calculate in system (|47|) the term 



(50) 



This prescription will come directly from the previous 
analysis and from the fact that matter keeps a frozen 
composition when it is perturbed. 



C. Motion of perturbed npe matter 

Let us first summarize the conclusions that come from 
the microphysical study overviewed in Section [llll 



- the unperturbed star is described by an effective 

barotropic EOS with the npe matter in beta equi- 
librium; 

- the perturbed (oscillating) star can not be described 

by the same EOS since matter is no longer in beta 
equilibrium; 

- the easiest but most general EOS for perturbed n-pe 

matter at low temperature depends on two param- 
eters: the baryonic density n?, and the proton frac- 
tion Xp; 

- for this nonbarotropic EOS and in the situation re- 

tained for the present study, the physical conditions 
in which the perturbations occur are that a given 
piece of npe matter will move with a frozen com- 
position that corresponds to its composition in the 
unperturbed situation. 

This last point needs now to be translated into an 
equation that will replace the beta equilibrium relation 
in the case of perturbed matter. But before we make 
it more explicit, we shall remind what turns out to be 
the situation when instantaneous beta equilibrium (or 
barotropicity) is assumed. This situation was supposed 
to occur in Paper I, and is still useful here since we have 
tried to modify our algorithms as little as possible. 

In this case, thermodynamics tells us that both the 
pressure P and the total energy density p are functions 
only of the baryonic density rif,. Hence, the term (|50|) 
can be written as the Eulcrian perturbation of an exact 
derivative (or the gradient of a scalar field), or equiva- 
lently as the gradient of the perturbation of a field that 
comes from an exact derivative, the relativistic enthalpy 



1 



/(rib) 



^P{ni,) ^S^^H^nb) = "^SHiub) . (51) 
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Using this new function, one can write the EE in 
a condensed and well-known form, which makes the 
system self-consistent and which was described in Paper 
I. 

When we can no longer use the same EOS for the per- 
turbed star and the background star, the previous list of 
conclusions are translated into equations in the follow- 
ing ways. First, as already mentioned, they imply that 
we need an EOS for the perturbed fluid that relates the 
perturbation of the pressure P and energy density p with 
(for instance) the perturbation of the baryonic density rif, 
and the perturbation of the proton fraction Xp 



6P{nb,Xp) 



dP 
dub 



6nh 



dP 

dXri 



SXu 



(52) 



and no longer only with the first one^ 



dP 

5P0{nb)^ — 
drib 



_ P 

Siib^ — 'fpSub, (55) 



as it was when instantaneous beta equilibrium and a 
barotropic EOS were assumed. Note that here we have 
introduced the adiabatic index for the matter in beta 
equilibrium, 7^. This index will be used in the following 
and is a priori not a constant (as it is in the polytropic 
case). Moreover, the function Xp'^{nb), which links the 
proton fraction with the baryonic density for npe matter 
in beta equilibrium, is of course the result (|ll|l obtained 
in Section ITU 

One of the main implications of an EOS dependent on 
two parameters is that it makes it impossible to define 
the relativistic enthalpy by the relation (|51|l . since the 
fraction is no longer an exact derivative. Hence, the 
term (|5U|) can no longer be regarded as a function of a 
single variable and the system has to be completed by 
an additional equation to be closed, an equation that de- 
scribes the fact that the matter has a frozen composition. 

This prescription means that the Lagrangian perturba- 
tion of the proton fraction is (as far as the EE are con- 



From the relativistic point of view, it is crucial to verify that the 
employed EOS dependent on these parameters is covariant. It 
can be shown to be trivial if we write the numbers ni, and Xp 
as Lorentz scalars defined with the baryon, neutron and proton 
4-currents. It gives 



where i is either b, n or p. Moreover, 
4-currents and 4-velocities 



(53) 

with the relations between 
(54) 



it is easy to see that assuming each of the velocity fields to be 
equal implies that the most general EOS can only depend on two 
Lorentz scalars, n;, and Xp for instance. 



cerned) or (equivalently) that its derivative along the ve- 
locity field vanishes. Addressing with A the Lagrangian 
perturbation of any quantity, we then have for the mo- 
tions that we are dealing with 



P 

AP = — -fpAnb. 
rib 



where 



7-F = 



d\n[P] 
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(56) 



(57) 



is the adiabatic index for a frozen composition. Further- 
more, the formal relation, which links Lagrangian and 
Eulerian perturbations by the Lie derivative with respect 
to the Lagrangian deplacement ^, 



enables us to arrive to the final expression 



6P = SP, 
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P A , 

— An;, (7_F 
rib 



(58) 



(59) 



where the first part of the perturbation is what we would 
get if beta equilibrium was instantaneous. A similar ex- 
pression can be found for the perturbation oi f = p -\- P 
with some "C" coefficients instead of the 7. Nevertheless, 
the relevant quantity at the linear order for matter close 
to beta equilibrium is 



P 

Cf -Cl3 = J ilF - 7/3) 



(60) 



The expression (|59ll and its equivalent for / are those 
we use in the EE in order to write the term (|50|) as the 
perturbation of the gradient of the enthalpy plus a con- 
tribution that exists only for nonbarotropic stars. In this 
way, we do not need to change our algorithms (see Pa- 
per I for more detail) for the solution of the EE and the 
anelastic equation (|49|l . whereas we only need to add, as a 
source term, the second part of these expressions. More- 
over, written like this, this source term depends only on 
the Lagrangian perturbation of the baryonic density Arib 
(as far as first order terms are concerned), since we have 

+ ^H^ (7i. -7/3) (1 - 7) : (61) 

in which Arib is obtained by solving the equation 

£uXp ^ 0, (62) 

where u is the 4- velocity. Indeed, when we linearize this 
equation and use the anelastic approximation, we get af- 
ter some algebra the Newtonian-like advection equation 
for Arib, 



{dt + nd^) Arib 



rib ■ 



(63) 
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where the usual Newtonian notations were used. This 
new equation enables the existence of modes related to 
the gradient of composition and whose restoring force is 
gravitation. These (7-modes for cold but heterogeneous 
NSs were first discussed by Reisenegger & Goldreich [s^ l 
and will be the subject of the next section. 



V. GRAVITY MODES 

A. General features of composition gravity modes 

Before dealing with the more complex situation of the 
modes of a stratified rotating relativistic star, it is worth 
to recall the situation for the modes of a nonrotating 
star, in order to better recognize in the following their 
scions. With the already mentioned assumptions, it 
is quite simple to describe. We consider only the npe 
core of NSs, without superfluidity, without magnetic 
field and without trying to follow modes of the crust. 
Hence due to the spherical symmetry, only polar fluid 
modes® can exist and the spectrum of axial modes 
is degenerate at zero frequency. Moreover, using the 
relativistic generalization of the anelastic approximation 
[cf. Eg. 1)491) . Paper I where it was introduced and also 
the discussion in the Appendix], the only hydrody- 
namical modes left in the nonrotating cold star are 
low-frequency modes whose restoring force is gravity, 
the so-called composition 3-modcs. Indeed, the main 
consequence of the anelastic approximation is to filter 
out pressure modes (see the Appendix). In addition, 
the relativistic Cowling approximation prevents us from 
having w-modes (some of which could have been axial), 
which are mainly oscillations of the spacetime itself. 

Gravity modes of NSs have already been the subject 
of several works, even if they were first thought to be 
degenerate_at zero frequency due to beta equilibrium 
(Thorne, |53|). Thus, gravity modes that exist because 
of the presence of temperature gradients in warm NSs 
were studied by Mc Dermott et al. 01 1 who introduced 
the relativistic Cowling approximation. They found that 
the thermal ^modes were concentrated into a thin layer 
close to the surface, as was the case in the analysis of 
Finn [s^ on gravity modes linked to discrete changes 
in composition of the crust. The first study considering 
Newtonian gravity modes associated with smooth 
composition gradients in the core of NSs is that of 
Reisenegger & Goldreich [s^- However, they neglect the 
effect of the strong interaction, which is inappropriate 
since nucleons are far from Fermi gases [see for instance 
Eq.l^nH- Later, Lai js^ studied the CFS instability of 



We remind that polar modes are those whose parity for a given 
spherical harmonic number / is (—1)' while axial modes have 
(— 1)'+^ for parity. 



composition gravity modes using EOSs that include the 
strong interaction (see H^), but in Newtonian rotating 
NSs. While Yoshida & Lee [H3| looked at the r-modes 
of relativistic stars with gravity modes, yet they mainly 
focussed on modes due to temperature gradients, and 
only dealt in a very approximative way with modes due 
to composition gradients. Going farther in tempera- 
ture, gravity modes due to discontinuity of density in 
warm nonrotating relativistic NSs were investigated by 
Miniutti et al. l58|, whereas even more recently gravity 
modes due to temperature and composition gradients 
in hot rotating relativistic NSs has been analysed by 
Ferrari et al. j59j, restricting themselves to polar modes. 
On the other hand, the situation of gravity modes in 
very cold (superfluid) NSs was first studied by Lee 
[60j, who realised that there was no gravity modes in 
models with two fluids. This was later confirmed by 
Andersson & Comer ^61] and Prix & Rieutord and 
is now considered as a crucial feature of superfluidity 
to prove its existence in NSs through the observation of 
gravitational wave signals. 

Nevertheless, as already explained, we shall deal here 
only with modes of not too cold neither too hot NSs, 
and the main characteristics of the relevant composition 
gravity modes resulting from the previously mentioned 
studies are 

(i) a weak coupling to the crust and the independence 

with respect to the transition radius between the 
core and the inner-crust (H^jH^)- Following this 
result, we decided that it was appropriate to ne- 
glect the coupling between the core gravity modes 
and the crust modes, just giving a boundary con- 
dition for the modes at the transition density. No- 
tice that in Paper I we found that using either free 
surface boundary condition or rigid crust boundary 
condition was not a key issue in the linear approach 
to inertial modes. Hence, in all the following, only 
modes with the boundary condition of a rigid crust 
(null radial velocity) at a given density are calcu- 
lated. This density is of course the chosen density 
for the cutoff of the PAL EOS (see Section IHD)! . 
i. e. it is either half the saturation density or almost 
two thirds of it. The influence of the cutoff on the 
mode shall be discussed in the next subsection; 

(ii) the growth time of their gravity driven instability, 

which is much longer than the time of viscous 
damping, makes the (7-modes irrelevant for the pro- 
duction of gravitational waves (at least for cold 
NSs, see ;59j). This comes from their weak cou- 
pling to the gravitational field. This consideration 
is pointless for a nonrotating star, but it explains 
why in the next section we shall only deal with 
the current quadrupole that makes m — 2 inertial 
modes unstable; 

(ill) their frequencies, which are roughly between 50 and 
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400 Hz for a M = IAMq star (|53,|53,|5a), de- 
crease for a slightly decreasing size of the core. 

The latter point turns out to be easy to understand 
from the usual dispersion law for (?-modes that can be 
reached in the Newtonian case using the WKB approach 
(see the Appendix) 



1(1 + 1) 



(64) 



where w is the frequency, I the quantum number of the 
decomposition in spherical harmonics, r the radial coor- 
dinate, k the radial component of the wave vector and 
M"^ the Brunt- Vaisala frequency defined as 



1 



|VP|^ (7, 



(65) 



Notice that these 7 are the Newtonian equivalent of 
the relativistic gamma coefhcients^*^ introduced in Sec- 
tion IIVI and also that in the Appendix this Newtonian 
formula is written with the sound velocities which are 
easily related to those coefficients. 

Looking at Equation H64(l . we see that for a given /, to 
make the radius of the star a little smaller (which means 
to increase the value of k without changing M) im- 
plies a decrease oiw in accordance with point (iii) above. 

Since N depends on the background star, the variation 
of the frequencies for changes of mass or of EOS are more 
subtle. To have an insight on this in our case, it is useful 
to look at the relativistic equivalent of Eq. (|65|l 



Pip + P) 



(66) 



where we restored explicitly the velocity of light c. The 
values of this frequency N for our models are depicted 
in Figure ^ whereas Fig|2] shows, as an indication, the 
gamma coefficients for the models with maximal central 
densities [A(4) and B(3)]. In Figd we see that the Brunt- 
Vaisala frequency (that is an upper-bound for the allowed 
frequencies) seems to be higher for softer EOSs [compar- 
ison between models A(4) and B(3) that have the same 
mass], and increases with the mass for a given stiffness. 
In the previous studies, this point did not seem very clear, 
since Reisenegger & Goldreich found that the frequency 
of gravity modes increases with mass (which is consistent 
with the increase of the Brunt- Vaisala frequency), but 
the result of Lai with fixed mass was that an apparently 
softer EOS (smaller radius) leads to a smaller Brunt- 
Vaisala frequency. Yet, it is worth pointing out that 



The difference between relativistic and Newtonian gamma coef- 
ficients is tliat tlie baryonic density nt, lias to be replaced with 
the mass density n in the Newtonian case. 
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FIG. 1: Brunt- Vaisala frequency for all the models used in 
this work up to the radius for which the density is half the sat- 
uration density. Note that since the frequency vanishes at the 
center of the star (null gradient of pressure) and monotonously 
increases with the radius, gravity modes will have an allowed 
propagation zone with a lower limit. 



Y coefficients for models A(4) and B(3) 
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FIG. 2: For the EOSs A and B, 7 coefficients as functions of 
the radius (in units of the cutoff radii) in the stars with max- 
imal central densities in order to scan the wider ranges as 
possible. Notice that the frozen coefficients are always higher 
than the coefficients at equilibrium, which is a necessary con- 
dition for stability. 



the properties of a NS depend not only on the compres- 
sion modulus, but also on the symmetry energy. Since 
Lai used two completely different EOSs, his results are 
probably not an indication of the influence of the stiff- 
ness on the Brunt- Vaisala frequency nor on the gravity 
mode frequency. This issue will be addressed in the next 
section. 
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Spectrum for the Lagrangian perturbation of density 
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FIG. 3: Fourier power spectrum of the time evolution of a 
Lagrangian perturbation of density on the equator for m — 2 
initial data. The total physical duration is 1.5 s. See also 
Figure0]for a logarithmic representation. 



B. Test of the code and gravity modes 

The first test that we did was to take one of the back- 
ground stars presented in Section Hi Dl [model A(21)] and 
perturb it with a small m — 2 Lagrangian perturbation 
of density 

Am = (l - r^) sin(i9)2 cos(z9) 

X [cos(2 0) -f sin(2 (j))] . (67) 

The result and the stability of the code are illustrated 
in Figures |3| and ^ that depict the Fourier spectra of 
the evolution. Figure ^ is in logarithmic scale as will 
be every similar figure in the following: in this way 
more modes are visible. With this model, the higher 
frequency gravity mode has a frequency around 150 Hz. 
This order of magnitude is in agreement with previous 
results (IE^jISs'I), but we shall discuss their values and 
how they depend on the background star in the following. 

In Figure 01 we also put another calculation done with 
the same background star, but changing the initial data. 
Instead of using initial data whose angular decomposition 
corresponds to one single azimuthal number (m = 2), we 
tried initial data with m = and m = 2 mixed: 



Anb = 



I - r2) sin(i9)2 cos(i9) 
X [cos(2 0) + sin(2 0) - 1] 



(68) 



It can be shown that while the first data were only 
composed of the (/ = 3, to = ±2) associated Legendre 
polynomials, the second data additionally contain 
{I = 1,TO — 0). Hence, this figure shows that, as 
expected in a spherical star (nonrotating and without 
anisotropic physics), the modes that can be seen do not 



depend on the azimuthal numbers present in the angular 
decomposition. Exception done of the fact that for a 
given value of m only modes with / > \m\ can exist. 
This is the reason why for (m = 0) an additional mode 
with the frequency ^ 88.8 (Hz) (plus others with very 
small frequencies) appears. 

Since our main goal is to study inertial modes with 
m — 2 (see Section IVI|I that are the most interesting 
from the gravitational waves point of view, we shall re- 
strict ourselves only to to = 2 initial data as soon as we 
shall deal with rotating stars. But for the time being, we 
shall keep the mixed initial data in the study of gravity 
modes without rotation. The next point that we shall 
discuss is the influence of the cutoff density. Following 
the already mentioned works on gravity modes in NSs, 
we expect that this cutoff should hardly have an influ- 
ence on the spectrum, the main influence being linked to 
the change in the size of the core. This point is clarified 
by Figure |S1 that shows, for the same background star 
as before, spectra corresponding to the time evolution of 
the Lagrangian perturbation of (baryonic) density, but 
also of the radial and components of the (Eulerian 
perturbation of) velocity for two cutoff values: 0.08 or 
0.1 fm~'^. Note that these spectra are reached with 
evolutions done at different places inside the star. With 
Figure El we are able to verify that to cut at a higher 
density (which means to make the star smaller) indeed 
decreases the frequency, but the change is quite small: 
some 3% in the worst case. Hence, we shall now restrict 
ourselves only to calculations done with a cutoff value 
of 0.1 fm~^. Moreover, we also see in this figure that 
there is a very good agreement between spectra reached 
from time evolutions of density or of a component of the 
velocity. In addition, the fact that not all gravity modes 
have a radial velocity (gravity modes can be either polar 
or axial) is illustrated, mainly with low-frequency modes. 

The next step in our calculation was to look for the 
influence on gravity modes of the physical parameters 
that describe the background star. This is depicted in 
Figure |^ where we draw the frequency of the gravity 
modes for all models presented in Table Thus, the 
three first graphs depict spectra for a relatively stiff mat- 
ter (compression modulus of 180 MeV) with increasing 
central density (from top to bottom) whereas the last 
two curves are for an EOS which is stiffer: compression 
modulus of 240 MeV. The first comment to make is 
that to increase the mass (the central density) for a 
given EOS shifts the spectra to higher frequencies. This 
result is in agreement with the increase of Brunt- Vaisala 
frequency found in Figure ^ 

Next, comparison between the third and the fifth 
curves also gives a dependence in agreement with Fig- 
ure ^ These two curves indeed correspond to stars with 
the same masses, but with different stiffness. Hence, we 
verify that stiffer EOSs [model B(3)] lead to lower fre- 



15 



Spectra for the Lagrangian perturbation of density 
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FIG. 4: Same Fourier power spectrum as in Figure |21 but in 
logarithmic scale. In addition is pictured the spectrum for 
another calculation that corresponds to a mixture of m = 
and m = 2 as initial data. The agreement concerning the 
values of the frequencies of the modes is very good and a 
I = 1 mode, which only exists for m = 0, is the main difference 
between the two spectra. 



quencies, even if they imply smaller stars. The radius 
rather than the mass seems to be the dominant param- 
eter to determine the value of the highest gravity mode. 
This is made more evident when we look at the second 
and the third curves [models A(26) and A(4)] and at 
the fourth and the fifth curves [models B(24) and B(3)]. 
Indeed, the difference between the mass of A(26) and 
A(4) is about 0.25 solar masses for a difference in the 
frequency of the highest gravity mode of 17 Hz (160 to 
177 Hz), while for B(24) and B(3), the difference between 
the masses is around 0.2 solar masses with only 5 Hz be- 
tween the frequencies (126 to 131 Hz). However, as we 
can see in Tabled the difference between the radii is of 
0.4 km for model A and only 0.03 km for model B. In 
fact, this result only illustrates the same as in Figure ^ 
matter added to a star with a stiff EOS hardly changes 
the radius. This is the reason why the Brunt- Vaisala 
frequency was changed between models A(26) and A(4), 
but was hardly modified between B(24) and B(3). Yet, 
as we shall see now, when the star is rotating it is more 
tricky to identify modes just by looking at a spectrum. 



VI. INERTIAL MODES AND COMPOSITION 
GRAVITY MODES IN ROTATING NEUTRON 
STARS 

A. Modes in rotating stars 

As explained in the Introduction, modes of rotating 
relativistic stars have the peculiarity that some of them 
can be driven unstable by their coupling to gravitational 



waves. In fact, this CPS C^,^) instability IS generic, 
and only insufficiently high angular velocity of a star 
or physical phenomena like viscosity can prevent some 
oscillations of being unstable. Furthermore, Friedman 
& Schutz proved that the instability appears when 
an initially retrograde mode becomes prograde (as seen 
by an inertial observer). To explain how we should be 
able to "detect" this occurrence in the present work (in 
which we only look at spectra with positive frequencies), 
it is worth to start with a brief summary of the main 
differences between modes of rotating and nonrotating 
stars. 

The first obvious effect of introducing rotation in the 
Euler equation for the oscillations of a star is to allow 
the appearance of modes restored by the Coriolis force 
(the so-called inertial modes). But rotation also breaks 
the degeneracy at zero frequency of purely axial modes. 
This leads to the possible existence of purely axial iner- 
tial modes, the so-called r-modes, while it also possibly 
couples together axial and polar modes. Another effect 
of rotation, on the modes that already exist in the non- 
rotating star, is the splitting of their frequencies. Thus, 
frequencies measured in the inertial Wi and rotating Wr 
frames are linked by the relation 

Wi = Wr — mfl , (69) 

where m is the azimuthal number and il the angular 
velocity of the star (or of the rotating frame). Since 
we now restrict ourselves to modes with |m| = 2 (the 
most interesting for GW), this splitting should appear 
in our spectra calculated in the inertial frame as the 
replacement of any mode (with frequency wq in the 
nonrotating case) by a pair of modes with frequencies 
Wr + 2Q and Wr — 217, with Wr wq. Notice that 
the small difference between Wr and wq is linked to the 
influence of rotation on the structure of the star and on 
the modes frequencies, a difference which should remain 
quite small since we keep here only terms linear in f2 
(slow rotation approximation). 

For such a pair of frequencies, the mode whose positive 
frequency decreases with increasing f2 is obviously the 
retrograde mode. Hence, in Section IVl CI where we shall 
have gravity modes in rotating stars, a possibly unstable 
gravity mode will be detected when in the spectra (for 
increasing f2) its frequency have reached zero and then 
started to grow again. But as mentioned in Sections 
nil Dl and gravity modes of cold NSs are not good 
candidates for the emission of gravitational waves due 
to the CFS instability. They are weakly coupled to the 
gravitational field and we shall then not pay too much 
attention to their potential instability. 

The situation is quite different for inertial modes. In- 
deed, restored (in the Newtonian limit) by the Coriolis 
force, they only exist for rotating stars and have frequen- 
cies proportional to the star's angular velocity fl. Thus, 
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FIG. 5: Fourier power spectra of Lagrangian density perturbation, radial and j? components of velocity for a star with two 
different values of the density at which is done the cutoff. As expected for gravity modes, the frequencies decrease when the 
size of the star is made smaller by increasing the value of the cutoff density. Yet, the difference is quite small. 



the Newtonian frequency of a purely axial inertial mode 
is, in the rotating frame, 



Wr 



2mn 



(70) 



which gives the always negative frequency (which means 
an always prograde mode) in the inertial frame 



i{i + i) 



2n 



1(1 + 1) 



(71) 

This is the reason why, as mentioned in the Intro- 
duction, some inertial modes are unstable whatever the 
angular velocity of the background star is. Moreover, 
it can be shown (I3|, [5|) that depending on the 
assumed motion (whether or not the perturbations are 
adiabatic), the properties of inertial modes change. 
Thus, a Newtonian barotropic star has a spectrum in 
which only axial modes (r-modes) subsist for spherical 
harmonics that satisfy / = m, whereas nonbarotropic 
stars admit axial modes for any combination of I and to. 

In the relativistic case, the situation is a bit different, 
starting with the fact that no pure axial modes persist 



for barotropic stars as proven by Lockitch et al. [63 (see 
also Paper I). For nonbarotropic relativistic stars, the 
problem was quite controversial during some time and it 
is still not so clear. Indeed, in the relativistic framework, 
purely axial modes of nonbarotropic slowly rotating stars 
are expected to satisfy the Kojima's master equation [6^ . 
This equation contains some perturbations of the metric, 
but it was shown 5] that the main features of r-modes do 
not change with the relativistic Cowling approximation. 
Moreover, with this approximation, due to the existence 
of the frame-dragging effect, this equation easily leads to 
the relativistic equivalent of Eq. H71|l 



Wi 



n 



2m 
1(1 + 1) 



in which 



r^n - N. 



if J 



(72) 



(73) 



where is the third component of the shift vector 
defined in the metric H42|) . Since w is a, function of r (and 
only of r in the slow rotation approximation), Ea. H72|) 
implies a "continuous spectrum" (see j65| for a proper 
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Spectra for gravity modes for different EOSs and masses 
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FIG. 6; Fourier power spectra of Lagrangian density perturbation for stars with different masses and EOSs. All the models 
used in this study are presented here. 



demonstration). This made unclear for some time the 
existence of discrete r-modes for certain relativistic stars 
(for more detail see the review by Kokkotas & Ruoff 
Q)- For instance, Ruoff and Kokkotas (see references 
within [3) found that only for a very restricted range of 
polytropic stars was possible to encounter purely axial 
inertial modes. Thus, the larger is the polytropic index 
of a star, the lower its maximal compactness can be for 
it to admit discrete r-modes in its spectrum. 

But it was also demonstrated 'G^I that Kojima's master 
equation becomes a singular eigenvalue problem if the 
frequency in the rotating frame, Wr^wt + m 51, verifies 



2mw{0) 



< Wr < 



2mzu{R) 
l{l + l) 



(74) 



where is the null radius (center of the star) and R the 
radius at the surface. For this reason, it was claimed by 
Lockitch & Andersson 66] that some "boundary layer 
like" approach could be the way to properly solve this 
problem. Later Ruoff et al. realized |63 that to include 
a coupling between I terms of axial modes with / ± 1 
terms of polar modes (and respectively) could imply 
the existence of modes with discrete frequencies, some 



of them possibly hidden in the continuous part of the 
spectrum. However, they also found that this continuous 
spectrum was hugely influenced by the number of / they 
coupled together. 

Thus, the oscillation spectrum of a nonbarotropic ro- 
tating relativistic neutron star is still an open problem 
that we shall deal with in Section FVI CI in which modes 
of the background NSs described in Table will be stud- 
ied taking into account the frozen composition. Hence, in 
order to check more easily if the frequencies that will be 
displayed are inside or outside of the range (|74|l . we sum- 
marize, in Table ITTl their limits for (m = 2 , ^ = 2 or 3) 
and for all our models (c/. Table Pi. But since these 
limits are defined in the rotating frame while we work in 
the inertial frame, we put in this table not directly w but 
defined as 



2 TO n7 



(75) 



Finally, to better recognize the influence of the assump- 
tion of frozen composition, we shall first focus on inertial 
modes of those stars with the hypothesis of barotropic 
EOSs. 



18 



TABLE II: Dimensionless limits of the ranges of continuous spectra of some purely axial inertial modes as seen in the inertial 
frame for all models presented in Tabled] In order to use these limits for arbitrary angular velocity Q, the dimensionless 
is introduced. 



EOS 


Model A(21) 


Model A(26) 


Model A (4) 


Model B(24) 


Model B(3) 


nsurf (fm"^) 


0.08 0.1 


0.08 0.1 


0.08 0.1 


0.08 0.1 


0.08 0.1 


tSj{0)32 


-1.602 -1.600 
-1.801 -1.800 


-1.643 -1.642 
-1.821 -1.821 


-1.731 -1.733 
-1.865 -1.866 


-1.631 -1.631 
-1.816 -1.816 


-1.681 -1.675 
-1.841 -1.837 


^{R)22 
tSj{R)32 


-1.422 -1.423 
-1.711 -1.711 


-1.436 -1.438 
-1.718 -1.719 


-1.468 -1.473 
-1.734 -1.736 


-1.436 -1.439 
-1.718 -1.719 


-1.459 -1.455 
-1.729 -1.727 



B. Modes of relativistic rotating barotropic stars 

The first thing that we shall look at is the influence 
of our cutoff density on inertial modes. To explore 
this, we did several time evolutions with the Newtonian 
I — m — 2 r-mode as initial data (see Paper I). This 
purely axial mode was evolved in all backgrounds with 
different angular velocities of the star, and a typical 
result is illustrated by FigH This figure depicts the 
spectra of the d component of the velocity for a star 
of model A(21) with an angular velocity of 50 rad.s"^. 
Three different values for the cutoff density are used: 
0.053, 0.08 or 0.1 fm~'^, that is to say more or less 1/3, 
1/2 or 2/3 of the saturation density rip. Moreover, in 
order to test the possible appearance of a continuous 
spectrum, we draw the spectra calculated at 2 different 
positions within the star. Thus, the first obvious result 
is that for this model at that precise angular velocity, 
we see a discrete spectrum, with the Z = m = 2 inertial 
mode (the higher peak in the spectrum) that seems to 
have a frequency around 73 Hz. This value corresponds 
to a ratio between the relativistic and Newtonian 
frequencies in the rotating frame which is around 0.82, 
in quite good agreement with the results of Lockitch 
et al. |68j | for that compactness. In addition, various 
inertial modes can be seen that can be identified as the 
relativistic counterparts of the inertial modes found by 
Lockitch & Friedman 69], but since it was shown |63| 
that none of these modes are relevant for GW with the 
CFS instability, we shall not describe them in detail. 

Moreover, as can be seen, the influence of the cutoff 
density on the spectrum of inertial modes is even 
smaller than for gravity modes. This can probably 
be understood from the fact that starting with purely 
axial initial data we reach a radial velocity smaller 
than in the case of the 3-modes. This last statement 
is supported by the comparison between Figures 
and |S| Both of these figures depicts spectra for time 
evolutions of m = 2 modes in a star of model A(21). 
However, for the first one, the spectra of the and 
r components of the velocity correspond to gravity 
modes in a nonrotating star [initial data identical to 
the Ea. (|57|l and still without velocity], while for the 
second they correspond to the Newtonian inertial mode 
in a barotropic A(21) star with 17 = 50 rad.s"^. Notice 



that for both calculations the cutoff density is 0.1 fm^^. 
What can be seen is that the orders of magnitude of the 
1? component and of the radial component for the main 
peak are almost the same for gravity modes (Fig[SJ), 
whereas in the case of inertial modes (Fig|HJ) the radial 
velocity is more than one order of magnitude smaller 
than the •§ component. This feature can of course also 
be explained by the fact that the spectra of relativistic 
inertial modes mainly depend on the compactness of the 
star as we shall further verify. Notice in addition that, 
as expected, this inertial mode is not a purely axial mode. 

The last curve of Fig|Sl is the spectrum of one of the 
two independent components of the current quadrupole 
tensor (see Paper I for the exact definition). Since this 
tensor dominates the post-Newtonian reaction force , 
we shall use it as an indicator of the fastest growing 
modes, for those that verify the CFS criterion. This 
indicator tells us that the already mentioned inertial 
mode actually seems to be the most unstable one. 
However, even if fewer modes can be seen here than in 
the spectra of the velocity's components, several other 
modes appear. Among them, a peak around 47 Hz but 
also a quite hidden one around 81 Hz that is probably 
the axial-led I — 4 mode whose stability was discussed 
by [68j ]: it has indeed a frequency in the rotating frame 
that is some 70 % of the frequency of the I = m = 2 
mode, while in the Newtonian case, Lockitch et al. 
[68| found a ratio of 69%. Notice that this argument 
relies on the assumption of "almost similarity" of the 
spectra calculated in the rotating frame for different 
stars. Formulated in another way, it just means that 
the spectra in the rotating frames mainly depend on 
one global physical parameter, the compactness. This 
assumption will be tested in the following. Anyway, this 
mode was shown not to be relevant for GW, which is 
supported by its very low presence (compared with the 
main peak) in the Sij spectrum. It was more noticeable 
in the velocity's spectra. 

Now that we have analysed the effect of the cutoff 
density, we shall be able to investigate the inertial modes 
for all our models. In Fig|3 one can see the same kind 
of spectra for the •& component of velocity using all our 
background models and changing the angular velocity 
from 10 to 80 rad.s""'^, while in Fig^| only the curves 
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Spectra for component of velocity 



Q. = 50 rad.s , model A(21) with various cutoffs 
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FIG. 7: Fourier power spectra of the i9 component of velocity for inertial modes in a barotropic star with various values of 
the cutoff density. The spectra are calculated at 2 different positions inside the star illustrating the fact that they are discrete 
modes. Numerous inertial modes can be identified. 



Spectra for inertial modes in a barotropic star 

n = 50 rad.s"' with model A(21) 
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for 17 = 80 rad.s" are depicted in order to better see 
the differences between different stars. Thus, the first 
comment is that the influence of the background star on 
the spectra is quite small. The higher resolution makes 
that this is better seen in Fig llOl but already in FigElone 
can observe a very small dispersion of the frequencies. 
Moreover looking carefully, one can discern a cross inside 
the highest square, i.e. the square denoting the I = 2 
inertial mode for model B(24). Thus, models A(26) and 
B(24) seem to have the same frequency for this mode, 
which supports the fact that more than the EOS, the 
compactness of the star plays a key-role in the difference 
between Newtonian and relativistic frequencies. Indeed, 
as can be seen in Tab. those stars have almost exactly 
the same compactness: 0.189 and 0.191 respectively. 



FIG. 8: Fourier power spectra of the i? and r components 
of velocity and also of one of the 2 independent components 
of the current quadrupole tensor Sij for inertial modes in a 
barotropic star with cutoff density ~ 0.1 fm~^. In addition 
to the I = m = 2 inertial modes, others can be seen. 



Next, to test the dependence with respect to ft of the 
frequencies, we shall divide all of them by the angular 
velocity of the star. Since it is pointless for gravity 
modes, we shall restrict this analysis to the I = 2 iner- 
tial mode, whose dimensionless frequency is illustrated 
versus the angular velocity for all compactness in Fig llll 
Nevertheless, before doing more comments about what 
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FIG. 9: Fourier power spectra of the •& component of velocity for m = 2 inertial modes using all our background models with 
different angular velocity. The star is assumed to be barotropic with a cutoff density at 0.1 fm~"^. 



turns out from this figure, it is worth to stress its 
limitation. Indeed, the version of the code used in this 
study works with a given time step of 0.01 ms, and the 
total duration of the evolutions is always of the order 
of Is, whilst the angular velocity changes. Hence, the 
error bars on frequencies expressed in Hz are of about 
1 Hz, but for frequencies normalised with the angular 
velocity fi, they are equal to lrad.s~^/fl. Thus, for 
ri = 25 rad.s"^, this is 4 % and this is less than 2 % for 
> 50 rad.s""'^. All of this explains why we did not 
put the data for = 10 rad.s^^ that were completely 
meaningless. 

On this other hand, higher angular velocities of the 
star imply higher frequencies that give better resolution, 
which support the compactness dependence of the fre- 
quency as shown by Fig llll However, for all angular 
velocities, the frequencies of the I — m = 2 mode for 
models with compactness 0.189 and 0.191 are the samc^^. 
Yet, our main purpose was not to study inertial modes 



of barotropic stars, but to try to have an idea of the 
influence of the microphysical conditions, which implies 
nonbarotropicity, on the spectra. This is the reason why 
instead of working with time step and duration scaled 
by the inverse of O (as we did in Paper I), we use here 
only "physical scales" that we shall also keep in the next 
section for the study of the modes of stratified rotating 
NSs. 



C. Modes of relativistic rotating nonbarotropic 
stars 

1. Time evolution of a density perturbation 

The first mentioned effect of rotation was the splitting 
of noninertial modes. Our investigation of the pulsations 
of rotating stratified relativistic NSs will begin with 
the illustration of this phenomenon within the study 



We verified that for higher frequencies inertial modes, the higher 
resolution enables us to make a distinction between them. We 



shall not discuss this since those modes are not relevant for grav- 
itational waves emission. 
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FIG. 10; Fourier power spectra of the 1} component of velocity for m = 2 inertial modes using all our background models for 
57 = 80 rad.s"^. The star is assumed to be barotropic with a cutoff density at 0.1 fm~^. 



of density perturbations in rotating stratified NSs. To 
simplify, we shall first restrict ourselves to the model 
A(21) with various angular velocities, but keeping the 
cutoff density at 0.1 fm^"^. In fact, since we have verified 
in the previous sections its low influence, we shall use in 
the following only that value of the cutoff density. 

FigEl depicts the spectra of the Lagrangian pertur- 
bation of density. The initial data consist of the same 
m = 2 perturbation of density [see Ea. Hfci7|) ] without any 
velocity as in Section IV Bl As expected, the splitting of 
the (;-modes can be seen, but in addition, we can observe 
the emergence of some probable inertial modes: see the 
peaks around 120 Hz and 150 Hz for SI = 80 rad.s"^. 
The main reasons why these modes can be thought to be 
inertial modes are that both their frequencies and their 
relative importance in the spectra grow with the angular 
velocity. Yet, to better describe them, it is worth to have 
a look at the spectra of the velocity in Figs^l and 1141 
where we show spectra of the and radial components 
of the velocity, respectively. 

First of all, Fig^l directly shows that not all modes 
were visible in the spectra of the density. Thus, in 



addition to the splitted gravity modes and to the already 
mentioned candidates inertial modes, there are many 
others that do not correspond to peaks in the density 
spectra. We shall now see how we are lead to classify 
them as axial-led inertial modes, following the classifica- 
tion of Lockitch et al. 16^. Obviously, they are inertial 
modes since their frequencies increase with the angular 
velocity. But our view of classifying them as axial-led 
is supported by Fig ll4l which shows that the radial 
velocity for those modes is very small comparatively to 
the radial velocity of others. Moreover, one of them even 
seems not to have any radial counterpart (at this scale 
of precision): the mode with a frequency of 116 Hz for 
17 = 80 rad.s^^, that is to say the mode with the lowest 
frequency among those axial-led inertial modes. 

Interestingly, Fig^J which depicts the spectra for the 
current quadrupole tensor, shows that this mode is the 
mode with the largest contribution to this tensor. Notice 
anyway that its amplitude is very low (logarithmic 
scale) due to the initial data. Nevertheless, for all the 
preceding reasons, we can now say that this is the 
I = m ~ 2 inertial mode, and that Table tells us 
that its frequency is inside the range of the continuous 
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Dimensionless frequency of the 1 = m = 2 inertial mode 

relativistic barotrops with different cotnpactness 
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FIG. 11: Dimensionless frequency versus angular velocity of 
the star for I = m = 2 inertial mode in barotropic stars of 
various compactness. Notice that due to the fact that evolu- 
tions lasted 1 s, the error bar is equal to 1 rad.s~^/f2. 



the existence of a "chemical" part of the physical energy, 
while the second is linked to fluxes through the surface. 
To calibrate the violation of the conservation due to 
fluxes, we did time evolutions of this kinetic energy in 
barotropic stars with cutoff densities 0.053, 0.08 or 0.1 
baryon.fm^'^. We found that for these values, the higher 
was the cutoff, the better was the conservation, which 
means that our usual cutoff at 0.1 baryon.fm^^ is the 
best of them. Notice that in all cases we saw temporary 
violations of the conservation, but the temporal mean 
value of the energy was always conserved. Moreover, 
even in the worst cases, the violation were never higher 
than 0.6%. On the other hand, testing the influence of 
the "chemical energy" , we verified that for nonbarotropic 
EOSs with a cutoff at 0.1 baryon.fm"'^, the mean value 
was no longer conserved, but the violation was only of 
the order of some 3 percents. For all these reasons, we 
shall work in the following with this (kinetic) energy as 
a not too bad indicator of the energetic behavior of the 
oscillations, at least as far as the coupling between polar 
and axial parts is concerned. 



spectrum. Indeed, its dimensionless frequency in the 
inertial frame is 1.45 whereas for model A(21) the range 
is [1.42,1.6]. Moreover, with the already displayed 
spectra of velocity, the spectra of a global quantity (the 
current quadrupole), plus various spectra of velocity 
calculated in several positions within the star^^, we 
verified that this is a discrete eigenvalue. 

Finally, before we try to investigate the I = m = 2 
inertial mode in more detail by changing the initial data, 
let us introduce some axial and polar "kinetic energies" 
that shall help us to study the coupling between polar 
and axial parts of the modes. What we shall call "kinetic 
energy" in the following was already introduced in Paper 
I and is defined as 

E =1 [ d^YiiW^ , (76) 

2 Jstar 

where n^ni, N[r] a[r] [see is the three 

dimensional vector field (the Eulerian perturbation 
of velocity) defined in Ea. H46|l and d'^V is the fiat 
elementary volume. 

With the anelastic approximation, this quantity can 
be shown to be conserved for a barotropic EOS and a 
vanishing density at the surface, conditions that were 
verified in Paper I. But for nonbarotropic EOSs with 
a cutoff density, two ways to violate this conservation 
exist. The first one, linked to the EOS, corresponds to 



We do not show the corresponding graphs here, but they are very 
similar to what was done in Fig|7|for the barotropic case. 



Indeed, we remind that our algorithm for the solution 
of Euler's equations (see Paper I) relies on the use 
of the radial velocity and two scalar potentials whose 
angular divergence and curl enable us to recover the 
usual spherical components. By using these variables, 
the separation of the polar and axial parts of the velocity 
field is done very easily, which naturally leads to the 
polar and axial kinetic energies. Thus, FigE| that 
displays the time evolution of the ratio between the 
polar kinetic energy and the total kinetic energy give a 
first insight on the energy flows between the polar and 
axial parts of the modes. 



In this figure, time evolutions of the ratio are shown, 
whose duration is only 200 ms in order to make visible 
the "quasistationnarity" that is reached after a while. As 
can be seen, starting with a purely polar perturbation, 
the proportion of the energy stored in polar modes is in- 
deed decreasing down to a fraction that is between 80 and 
90%, depending on the angular velocity of the star: the 
faster it rotates, the more energy fiows to axial modes. 
Notice anyway that since the initial data are a perturba- 
tion of density without any velocity, at the very begin- 
ning, there is, rigorously speaking, no kinetic energy. But 
in fact, from the very first time steps the perturbation of 
density implies a velocity which is initially purely po- 
lar. Moreover, we encountered quite similar results when 
we did some tests with as initial data a purely polar ve- 
locity field without any density perturbation. Only the 
ratio was slightly changed, whilst the time scales were 
identical. Nevertheless, the relevant issue concerning en- 
ergy fiows will be how barotropicity affects this evolution. 
This will be one of the topics of the next section. 
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Model A(21), m = 2, cutoff = 0.1 fm" 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s' 
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FIG. 12; Spectra of the Lagrangian perturbation of density for several angular velocities. The star is nonbarotropic and of 
model A(21), while the initial data are a m = 2 density perturbation. See also Figs. 1131 114l and ll5l 



2. Time evolution of the I = m = 2 mertial mode 

We shall now study the time evolution of oscillating 
out-of-equilibrium NSs with different initial data: the 
1 — 171—2 inertial mode instead of a density perturba- 
tion. In this way, the differences with the barotropic case 
shall be easier to notice, while we should also be able 
to see how much our previous results for nonbarotropic 
stars depend on the initial data. To start with, FigllTI 
depicts the spectra of the Lagrangian perturbation of 
density for this case, with the results already presented 
in Fig ll2l for comparison. The first thing to comment 
on is that even for low angular velocity, gravity modes 
appear, with, as expected, an amplitude that grows with 
il. But it shall also be pointed out that the I = 2 and 
I — 4 inertial modes are visible, even if we consider the 
density perturbation spectra, and that their amplitudes 
stay the same for all values of SI. But the coupling 
between axial and polar parts of the velocity field should 
be better investigated with the spectra of the velocity 
components. These are displayed in Figsll8l and 1191 

FigEl shows that the d component of the velocity is 
hardly affected by the fact that the star is a barotrop. 



The main difference between the spectra is indeed only 
that, for nonbarotropic stars, gravity modes appear, with 
amplitudes that become higher for faster rotating stars. 
Even the frequency of the I = m = 2 inertial mode is 
almost exactly the same. But on the other hand, the 
radial part of the velocity fFig ll9|l is completely changed 
and quite strongly depends on the angular velocity, 
which was predictable. Thus, for slowly rotating stars, 
gravity modes dominate and the radial velocity of the 
inertial mode is some two orders of magnitude lower 
than in the case of a barotropic star. Yet, for increasing 
angular velocity, the amplitude and the frequency of the 
radial part of the I — m — 2 inertial mode become 
more and more similar to what they are in a barotropic 
star. For SI = 80 rad.s" ^, there is less than one order 
of magnitude difference between the barotropic and 
the nonbarotropic cases for the I = 2 mode (whose 
frequency in the inertial frame is Wi ~ 116 Hz in this 
case), whereas there was two orders for S7 = 10 rad.s~^. 
On the other hand, it should be pointed out that 
also the amplitude of the ^-modes grows with SI 
(probably meaning that their "inertial" part is on the 
way to dominate buoyancy) while many of the inertial 
modes that existed for barotropic stars have disappeared. 
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Model A(21), m = 2, cutoff = 0.1 fm" 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s 
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FIG. 13: Spectra of the i9 component of velocity for several angular velocities. The star is nonbarotropic and of model A(21), 
while the initial data are a m = 2 density perturbation. See also Figs. I12lll4l and ll5l 



If we now look at FiglJOl that depicts the spectra for 
the current quadrupole tensor, we shall observe that 
from the point of view of GW emission, barotropicity 
seems not to have a major impact, at least in the 
spectrum. Indeed, the spectra of this tensor do not show 
remarkable differences, with only a small shift in the 
frequency of the probably most unstable mode that can 
be pointed out. Wc shall now investigate in more detail 
how this shift in the frequency depends on the EOS, 
after a brief comment on the time evolution of the ratio 
between polar and axial kinetic energies. 

Indeed, while all previous spectra were very compara- 
ble, the time evolutions of the polar and axial parts of 
the kinetic energy are quite different for barotropic and 
nonbarotropic stars. The latter situation is illustrated 
by Fig|2] (out-of-equilibrium case), and by Fig[221 
(barotropic case)^'^. Thus, one observes that, for a 



barotropic star, energy is fast flowing back and forth 
from the axial part to the polar part, leading to a kind 
of "stationary state" with a medium value of the ratio 
around 2%. On the other hand, for nonbarotropic stars, 
the coupling between polar and axial parts seems to 
depend much more on the angular velocity. Hence, for 
il — 10 or 25 rad.s""'^, even after 600 ms of evolution, 
the "final state" is still not reached and the ratio still 
grows. But the curves for higher angular velocities 
suggest that they should not get too high since the 
"stationary situations" obtained in those cases are 
characterized by ratios lower than for barotropic stars 
(always less than 1 %). 

The difference between the typical times is quite easy 
to understand. Indeed, energy flows should have for 
typical time scales the periods of the modes measured in 
the rotating frame. For nonbarotropic stars, whatever 
the angular velocity, some (^-modes exist with periods 
of some ms, which are then responsible for the quite 



Notice that evolutions in the barotropic case axe, as expected, 
analogous to the result presented for a relativistic polytrop in 



Paper I. 
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FIG. 14: Spectra of the radial component of velocity for several angular velocities. The star is nonbarotropic and of model 
A(21), while the initial data are a m = 2 density perturbation. See also Figs. 1121 113l and ll5l 



fast flows. However, in the barotropic case, only inertial 
modes can be observed. Since they have periods propor- 
tional to the inverse of the angular velocity, this is the 
explanation of the two features of the barotropic case: 
the self-similarity (for different values of O) of the time 
evolutions and the acceleration of the flows with fi. 

On the other hand, with some caution due to our 
numerous approximations (and also to our definition 
of "kinetic energy"), we could say that the result 
concerning amplitude of the ratios can be seen as 
an indication than stratification probably helps the 
unstable inertial mode (the I — m = 2 r-mode) to 
be driven to instability by GW since its coupling to 
polar modes (which are not driven to instability by the 
current quadrupole) is smaller. Yet, the answer to that 
question requires a proper treatment of the emission of 
gravitational waves without the Cowling approximation, 
while different viscosities for polar and axial modes 
should also complicate the problem. 



latter will be the easier: as already discussed, due to the 
error bars for the normalized frequencies, we are not re- 
ally able to see any difference. Moreover, since our study 
is linear in f2, we expect that the normalized frequency of 
the inertial modes should only depend on the compact- 
ness of the stars but not on f2. This view is supported 
in the barotropic case and in the nonbarotropic case by 
(respectively) both Figsllllandl23l Thus, the more com- 
pact a NS is, the higher will be its dimensionless fre- 
quency as measured by an inertial observer. However, 
the comparison between these two figures shows that tak- 
ing into account the fact that npe matter should keep a 
frozen composition, which implies a nonbarotropic EOS, 
leads to lower frequencies. Moreover, it is worth point- 
ing out that we verified that, for all of our models, the 
I = m = 2 inertial modes had discrete spectra^^, even 
if FiglJSl and Tab. ^1 demonstrate that all these modes 
were found within the range of the continuous spectra. 
Nevertheless, as already mentioned and as can be seen 



To conclude this work, we shall give an overview of how 
the frequency of the I — m = 2 inertial mode depends 
on compactness, barotropicity and angular velocity. The 



which was easier to see by looking at the current quadrupole 
tensor that is a global quantity. 
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Model A(21), m = 2, cutoff = 0.1 fm" 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s' 
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FIG. 15: Spectra of the current quadrupole tensor S for several angular velocities. The star is nonbarotropic and of model 
A(21), while the initial data are a m = 2 density perturbation. See also Figs. 1121 113l and ll4l 



in Fig ll9l they are modes with a radial velocity, which 
means that they are not purely axial. For very slowly 
rotating stars, this radial component is quite small, but 
it becomes as large as what it is for nonpurcly axial iner- 
tial modes of barotropic stars when the angular velocity 
increases. Finally, notice that even if models A(26) and 
B(24) have the same compactness and then the same fre- 
quency for the I = m — 2 inertial mode, due to the 
mixing with gravity modes that depends on the detail 
of the microphysics, their full spectra are quite differ- 
ent. This issue is illustrated by Figl^that depicts the 
spectra of the radial component of the velocity for those 
models with the out-of-equilibrium hypothesis. Hence, 
this probably means that more precise studies including 
more coupling between modes (nonlinear terms in or 
nonlinear corrections to the Euler equations) should en- 
able us to make some differences between stars of the 
same compactness, and then to start to probe the inner 
composition of NSs when gravitational waves emitted by 
their instabilities are detected. 



VII. CONCLUSION 

Inertial modes and gravity modes of slowly rotating 
relativistic stratified neutron stars have been investi- 
gated in the hnear regime. Our study is based on a 
spectral three dimensional evolutionary code which has, 
at the moment of writing, the restrictions of using the 
Cowling and the anelastic approximations. However, 
the effect of the latter approximation is mainly to 
kill high frequency acoustic modes that should not be 
strongly coupled to inertial and gravity modes for slow 
rotation. Additionally, the slow rotation approximation 
is reasonable for observed pulsars as, even for the fastest 
known pulsar, the pertinent dimensionless parameter 
is of the order of some percents. The validity of the 
Cowling approximation is more subtle since without 
perturbations of space-time we do not directly include 
gravitational waves in our calculations, whereas they 
are one of the main motivations for this work. Never- 
theless, previous studies have proven that the properties 
of inertial modes were not too much affected by its use . 



Another major feature of our work is that we tried to 
simulate as well as possible the microphysical conditions 
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Model A(21), m = 2, cutoff = 0.1 fm 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s" 
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FIG. 16: Time evolution ol the ratio between the polar kinetic energy and the total kinetic energy (as measured in the rotating 
frame) in the same star as in previous figures: a nonbarotropic star of model A(21). 



that should occur inside not too massive and not too 
cold neutron stars. Thus, by using a quite reahstic 
equation of state for nuclear matter valid even when 
the beta equilibrium is broken, we did the first three 
dimensional time evolution of inertial modes in stratified 
relativistic stars, taking into account the frozen com- 
position of every perturbed lump of npe matter. With 
this approach, we were able to show that the coupling 
between polar and axial modes quite strongly increases 
with the angular velocity. This happens in such a 
way that, in nonbarotropic stars, initially axial inertial 
modes develop a radial part that can be as large as the 
polar part of the general inertial modes of barotropic 
stars. Yet, several of the polar-lead inertial modes of 
barotropic stars do no longer exist in nonbarotropic 
stars. Moreover, we found that the characteristic time 
scale for the exchanges of energy between polar and 
axial modes is different for barotropic and nonbarotropic 
stars. Thus, the coupling of inertial modes with gravity 
modes in nonbarotropic stars seems to give birth to fast 
energy interchange between the polar and axial parts 
of the fluid motion. This phenomenon could have some 
implications on the way viscosity or gravitational waves 
emission act, but it would of course need further studies 



to draw robust conclusions. 

However, the comparison of inertial modes in stars 
with and without the assumption of frozen composition 
has already lead us to the conclusion that nonbarotrop- 
icity decreases the influence of compactness on the fre- 
quency of the r-mode, making in addition this frequency 
slightly lower in a nonbarotropic star. If further studies 
with different equations of state and composition support 
this result, it could mean that the measurement of the 
global parameters of a neutron star (compactness, mass 
and radius if absorptions lines are observed, binary sys- 
tem parameters, etc.) in conjunction with the observa- 
tion of gravitational waves emitted by some instabilities 
of the inertial modes would be very instructive about 
the inner structure of neutron stars. But such a Grail 
would probably need to first improve our understanding 
of the physics with wider explorations of the physical pa- 
rameters: centrifugal terms and deformation of the star, 
magnetic field, differential rotation, nonlinear coupling, 
etc. Among them, superfluidity is one of the key phe- 
nomena, since several studies proved that gravity modes 
are suppressed by its existence. As a consequence, even if 
superfluidity in old neutron stars is already strongly sup- 
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Model A(21), m = 2, cutoff = 0.1 fm" 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s' 
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FIG. 17: Spectra of the density perturbation for several angular velocities of the star for an initial I — m = 2 inertial mode 
evolving in a nonbarotropic star of model A(21). 



ported by both theory and observation (ghtches), the de- 
tailed analysis of the gravitational wave spectra of oscil- 
lating neutron stars could be the possibility to definitely 
demonstrate the superfluid nature of the neutron stars 
content by the lack of well understood features linked 
with gravity modes. 



APPENDIX A: DISPERSION RELATIONS 

In this Appendix, the usual basic way to derive the 
Newtonian dispersion relation for 5-modes is summa- 
rized, and we briefly discuss the effect of the divergence 
free and anelastic approximations. 
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The easiest way to reach the dispersion relation for 
gravity restored modes is based on the perturbative La- 
grangian approach. If we restrict ourselves to linearized 
equations, and write ^ the displacement vector, A the 
Lagrangian perturbations and 6 the Eulerian (linked by 
the symbolic relation A^S + C*'^), we have for the Euler 
equation 



Sn 
11 



and for the mass conservation equation 
An = -n^-i, 



(Al) 



(A2) 
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Model A(21), m = 2, cutoff = 0.1 fm" 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s' 
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FIG. 18; Spectra of the i? component of velocity for several angular velocities in a barotropic and nonbarotropic star of model 
A(21) with inertial mode for initial data. 



where P is the pressure, (j) the gravitational potential 
and • indicates the usual Euclidean scalar product. 

If we now consider only spherical nonrotating stars 
(only radial dependence for (j)) and neglect the Eulcrian 
perturbation of the gravitational field (Cowling approxi- 
mation) , the nonradial part of Eq.imi) and Ea. jX2|) give 
us 



An 



(A3) 



in which we assumed a decomposition in spherical har- 
monics of the scalar variables [I is the index of the spheri- 
cal harmonic Yf^li), (p)] and a harmonic time dependency 
of^: 

dU^~w'l (A4) 

Introducing now the equilibrium Cf,q and frozen com- 
position (adiabatic) cp speeds of sound, 
|VP|- "2 

and 



AP = cl An , 



with the Brunt- Vaisala frequency defined as 



(A5) 
(A6) 

(A7) 



we obtained the following system of equations 



^d.{rH.)+^(l-^')-^4 = (A8) 
r-' ^ ' nCp \ w J Cp 

}_dr5P + — -(M' ^W') = 0, (A9) 

n Cp n 

with the Lamb frequency defined as wq'^ = cj^ l{l + 1) / r'^ 
with g the local gravitational acceleration. 



The new variables 



with 

A= - 

enable to write the system as 



(AlO) 
(All) 

(A12) 
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Model A(21), m = 2, cutoff = 0.1 fm" 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s' 
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FIG. 19: Spectra of the radial component of velocity for several angular velocities in a barotropic and nonbarotropic star of 
model A(21) with inertial mode for initial data. 



drV 



r^e^^ {wll w'^ - l) / (ncp) 



(A13) 



r 



with 



X 
Y 



A WKB-like approximation, 

drV = ikV, 
gives the dispersion relation 



1 I w, 



w 



Since for neutron stars we have 

TV < — < cj^/c, 

CF 



(A14) 



(A15) 



(A16) 



(A17) 



Ea. HA16|l leads to the two asymptotic and separated 
branches for the p-modes 



2 _ CF- 



dkr)^ + (A18) 



and the g-modes 



w 



[krY + /(/ + 1) 



(A19) 



2. Divergence free and anelastic approximations 

The two approximations that we shaU now discuss dif- 
fer by the way they treat the continuity equation ljA2|l . 
The divergence free approximation consists, as indicated 
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Model A(21), m = 2, cutoff = 0.1 fm 

From top to bottom: Q. = 10, 25, 50, 70 and 80 rad.s' 




FIG. 20: Spectra of the current quadrupole tensor for several angular velocities in a barotropic and nonbarotropic star of 
model A(21) with inertial mode for initial data. 



by its name, in the replacement of it by the condition In the case of the anelastic approximation, some alge- 

bra shows that the system (|A13p is replaced with 

V-S^ = , (A20) 

while the anelastic approximation neglects the time vari- 
ation of the Eulerian perturbation of density, which gives 

^■(ntj = 0. (A21) 



drV 







V . 



(A22) 



which leads to the relation IjAlQjl as the unique solution: 
the p-modes have been "killed" by the anelastic approx- 
imation. 

On the other hand, it turns out that the same exercise 
with the divergence free approximation also gives the re- 
lation (jA19|) as the unique solution, but the vector on 



which the WKB approximation is done has to be defined 
as 

This small difference in the definition of the second 
component of the vector between the divergence free ap- 
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Model A(21), m = 2, cutoff = 0.1 fm 
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FIG. 21: Time evolution ol the ratio between the polar kinetic energy and the total kinetic energy in the same nonbarotropic 
star as in previous figures with inertial mode for initial data. 



proximation and the "exact solution" seems to indicate such a difference) is closer to the exact solution than the 
that the anelastic approximation (for which there is not divergence free approximation is. 
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Model A(21), m = 2, cutoff =0.1 fm (barotropic case) 
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FIG. 22: Time evolution of the ratio between the polar kinetic energy and the total kinetic energy in the same barotropic star 
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FIG. 23: Dimensionless frequency versus angular velocity of 
the star for I = m = 2 inertial mode in nonbarotropic stars 
of various compactness. Notice that due to the fact that evo- 
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FIG. 24: Spectra of the radial velocity for several angular velocities in nonbarotropic stars of model A(26) and B(24) with 
inertial mode for initial data. 
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